Development of Monte Carlo configuration interaction: 
second-order perturbation theory 
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Approximate natural orbitals are investigated as a way to improve a Monte Carlo configuration interaction 
(MCCI) calculation. We introduce a way to approximate the natural orbitals in MCCI and test these and 
approximate natural orbitals from MP2 and QCISD in MCCI calculations of single-point energies. The 
efficiency and accuracy of approximate natural orbitals in MCCI potential curve calculations for the double 
hydrogen dissociation of water, the dissociation of carbon monoxide and the dissociation of the nitrogen 
molecule are then considered in comparison with standard MCCI when using full configuration interaction 
as a benchmark. We also use the method to produce a potential curve for water in an aug-cc-pVTZ basis. A 
new way to quantify the accuracy of a potential curve is put forward that takes into account all of the points 
and that the curve can be shifted by a constant. We adapt a second-order perturbation scheme to work with 
MCCI (MCCIPT2) and improve the efficiency of the removal of duplicate states in the method. MCCIPT2 is 
tested in the calculation of a potential curve for the dissociation of nitrogen using both Slater determinants 
and configuration state functions. 



I. INTRODUCTION 

Monte Carlo configuration interaction (MCCI) 1,2 of- 
fers the prospect of capturing many of the aspects of the 
full configuration interaction (FCI) wavefunction but us- 
ing only a very small fraction of the configurations. The 
method repeatedly performs a configuration interaction 
calculation with a set of determinants which is enlarged 
by the addition of random single and double substitu- 
tions and reduced by the removal of states which have a 
coefficient in the wavefunction of magnitude less than a 
specified cut-off (c m i n ). In principle the method can be 
applied to ground and excited states of single-reference 
or multi-reference systems. 

Single-point energies have previously been calculated 
using MCCIji as have the bond dissociation energies of 
water and HF.— The errors of electronic excitation ener- 
gies for atoms computed using MCCI were found to be 
small when compared with experiment in Ref. |j. Elec- 
tronic excitation energies for small molecules have also 
been calculated using MCCI with errors of generally circa 
ten meV when compared with experiment yet using only 
a tiny fraction of the FCI space. 5 Potential curves have 
been calculated for a variety of small systems using MCCI 
where it was found that non-parallelity errors approach- 
ing chemical accuracy when compared with FCI results 
could be produced using only a very small percentage of 
the FCI space even when the system was multireference^ 
However the calculation of the curve for the very chal- 
lenging system of fifty hydrogens presented difficulties 
for the method. Multipole moments, a non-variational 
quantity, have also been demonstrated to be calculated 
satisfactorily by MCCI for ground and excited states us- 
ing only a very small fraction of the FCI spaced MCCI 
ionisation energies for atoms have been shown to com- 
pare favourably with FCIQMC and exact results, while 
electron affinities were more challenging with larger per- 
centage errors but the absolute difference was not so 



poorii Again the size of the MCCI space tended to be 
very small compared to the almost always computation- 
ally intractable size of the FCI space. 

In this work we consider improving a MCCI calculation 
by using approximate natural orbitals or second-order 
perturbation theory. Although the Hartree-Fock molec- 
ular orbitals give the lowest energy single Slater deter- 
minant they may not be the most efficient choice for a 
CI calculation. The natural orbitals (NOs) are the eigen- 
functions of the first-order reduced density matrix or one 
matrix^ which are considered to give better convergence 
than Hartree-Fock molecular orbitals. One possible ben- 
efit is that some natural orbitals may have eigenvalues 
(occupations) which are essentially zero so can be dis- 
carded and hence the size of the FCI space is reduced. 
For methods where a wavefunction is not easily available 
or defined, derivatives may be used to calculate the re- 
sponse or relaxed density matrix which can then be used 
to give an approximation to the natural orbitals. It has 
been demonstrated that the NOs are indeed optimal for a 
system of two electrons^ but it is not clear if they are al- 
ways the best choice for larger numbers of electrons and 
Ref. [JJJ suggests that split-localised orbitals may offer 
better convergence in larger systems. Approximate nat- 
ural orbitals have been investigated as a possibly more 
efficient alternative to variationally optimising orbitals in 
CASSCF in Ref. [O for a CASCI calculation. There it 
was found that potential curves, including the dissoci- 
ation of ethylene, produced using CASCI with natural 
orbitals tended to usually have a non-parallelity error of 
only a few kcal/mol compared with the CASSCF curve. 
The exception was the approximate natural orbitals from 
restricted MP2 which tended to perform poorly while 
the best results were achieved with CCSD approximate 
natural orbitals. Natural orbitals from CISD calcula- 
tion for the ground-state of higher spin states have also 
been used for excited state MRCI calculations in Ref. [T^- 
There excitation energies were found to have a differ- 



ence of around 0.1 eV when compared with results using 
CASSCF/MRCI. 

We investigate the ability of approximate natural 
orbitals to improve the efficiency of an MCCI cal- 
culation. We use approximate natural orbitals from 
Quadratic CI with single and double substitutions 
(QCISD) 13 or second-order M0ller-Plesset perturbation 
theory (MP2).— For a multireference system these meth- 
ods may perform poorly or even fail to give sensible re- 
sults with regards to the energy so we also consider ap- 
proximate natural orbitals from an MCCI calculation. As 
a fairer comparison than just a single energy calculation, 
we consider if MCCI with natural orbitals can offer im- 
provements in accuracy and calculation time compared 
with standard MCCI for potential curves of water, car- 
bon monoxide and the nitrogen molecule when using FCI 
as a benchmark. 

We saw in Ref. [g that potential curves for small sys- 
tems for which full configuration interaction results were 
available could generally be calculated to relatively high 
accuracy using MCCI. This was achieved with a very 
small fraction of the states and it is interesting to con- 
sider whether results can be improved by using the MCCI 
wavefunction as the starting point for a second-order per- 
turbation calculation. To this end we adapt a second- 
order multireference perturbation method^ to work with 
MCCI and improve the efficiency of the removal of dupli- 
cate states. This method estimates the energy contribu- 
tion from the neglected states in a MCCI wavefunction at 
the expense of the final energy not being variational nor 
being easily associated with a wavefunction. We test the 
assumption that the MCCI wavefunction will be a very 
good starting point for this perturbation so that MC- 
CIPT2 should be able to produce a more accurate po- 
tential curve than MCCI, when both are compared with 
FCI, by accounting for more of the neglected dynamic 
correlation from a MCCI calculation. 



II. METHODS 



A. MCCI 



In this work we also consider a version of MCCI with 
a modified behaviour for the removal/addition of states, 
a convergence criterion and we also use Slater determi- 
nants (DETs) instead of CSFs for some computations, 
including the calculation of the MCCI NOs. When us- 
ing DETs the MCCI wavefunction is not necessarily an 
eigenfunction of S 2 and more states may be required but 
the construction of the Hamiltonian matrix and first- 
order reduced density matrix is much simpler. We cal- 
culate the Hartree-Fock molecular orbital integrals using 
MOLPRO J£ In this work the initial MCCI wavefunction 
is the CSF or DET formed from the occupied restricted 
Hartree-Fock orbitals. 



B. Natural orbitals 

The first-order reduced density matrix or one matrix 
is defined as 



j(xa,x b ) =N ^*(x a ,x 2 ,---,xn) 

^(xb,X2, ■ ■ ■ ,x N )dx 2 ■ ■ -x*n (1) 

which can be written in terms of the M one-particle 
molecular orbitals as 



M M 
i=l j=l 



(2) 



We use MCCI with Slater determinants and a not too 
onerous number of iterations with the same c m i n as the 
full MCCI calculation to construct approximate natural 
orbitals. We create the one matrix in the one-particle 
representation using the following method, beginning 
with 7 = 0. We consider all the DETs forming the wave- 
function. DETs i and j in maximum coincidence only 
contribute to 7 if they either have no differences to give 
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(3) 



where m runs over all orbitals in the DET, or one differ- 
ence due to orbitals k and I which results in 



The algorithm^ for MCCI is that the current MCCI 
wavefunction (usually initially comprising the occupied 
Hartree-Fock orbitals) has configuration state functions 
(CSFs) consisting of random single and double substitu- 
tions added to it. These substitutions are definitely at- 
tempted in CSFs with a coefficient of magnitude greater 
than a certain value while other CSFs have a 50% chance 
of a substitution occurring. The Hamiltonian matrix and 
overlap matrix are then constructed and the new wave- 
function is found. Newly added states whose absolute 
value of coefficient is less than c m i n are discarded and 
the process continues. Every ten iterations all states are 
considered for removal, not just newly added ones. This 
also occurs on the second last step and no states are 
added or removed on the final iteration. 



7/cZ ~~ ► Ikl 



&z>Ci Cj 



(4) 



Here e p is the sign due to putting the Slater determinants 
in maximum coincidence. We average over spins and then 
diagonalise the one matrix to give the MCCI natural or- 
bitals. As we cannot be sure that a very small occupation 
is due to the approximation rather than being something 
that would occur in the FCI natural orbitals, then we 
include all the approximate natural orbitals. We recal- 
culate the one and two-electron integrals using these ap- 
proximate natural orbitals in MOLPRO 16 then use them 
in a longer MCCI calculation using either DETs or CSFs. 
We also consider the approximate natural orbitals from 
QCISD and MP2 calculations. QCISD can perhaps be 
thought of as a less complex approximation to coupled 



cluster singles and doubles (CCSD)^ In QCISD size 
consistency is introduced into a configuration interaction 
method but at the expense of the energy no longer being 
variational. MP2 uses the Hartree-Fock Hamiltonian as 
the zeroth-order approximation in a second-order pertur- 
bation to give an efficient way to account for some of the 
correlation. It is also size consistent but not variational. 
Using MP2 or QCISD the one matrix may be approxi- 
mated by the response one matrix to give approximate 
natural orbitals which we generate with MOLPRO J£ Al- 
though the natural orbitals we consider are approximate 
we shall refer to them as the natural orbitals from a cer- 
tain method, e.g., MCCI natural orbitals. 



III. SINGLE-POINT CALCULATION USING NATURAL 
ORBITALS 

We first consider carbon monoxide at its experimen- 
tal equilibrium geometry^ with a cc-pVDZ basis, c m i n = 
5 x 10~ 4 , and the two lowest energy MOs or two most 
occupied NOs frozen. We use 500 iterations of MCCI 
with Slater determinants, but we see in Fig. [T] that the 
calculations have essentially converged in much fewer it- 
erations. The MCCI natural orbitals are calculated using 
a fifty iteration MCCI run. 



lowest energy (—113.042 Hartree) was from using QCISD 
NOs followed by MP2 NOs while the highest was for MOs 
(—113.036 Hartree). In Table |l] we display the time and 
number of DETs required to reach this energy. The time 
for the initial creation of the one and two electron inte- 
grals is not included as it is the same for each method. 
For MCCI NOs the integrals need to be recalculated and 
the time cost of this is included as is the QCISD and MP2 
calculation time. We see that MP2 NOs took the least 
time while QCISD NOs used the fewest number of DETs. 
MCCI NOs were an improvement of the time and num- 
ber of DETs compared with MOs but performed less well 
than the other two types of approximate natural orbitals 
we considered. As a system moves away from equilibrium 
it may be that the most efficient NOs are not from the 
same method. 

TABLE I. Total time and number of DETs for CO to reach 
E < —113.036 Hartree on the step following when all states 
have been considered for removal. 



Orbitals 


DETs 


Time 


(seconds) 


MOs 


9,919 




248 


MCCI NOs 


9,019 




173 


MP2 NOs 


7,453 




115 


QCISD NOs 


7,272 
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FIG. 1. MCCI (c m i„ = 5x 10 -4 ) energy against number of 
iterations for CO at R = 2.1316 Bohr with a cc-pVDZ basis 
set and with two frozen core orbitals when using either MOs, 
MP2 NOs, QCISD NOs or MCCI NOs. 

There is substantially faster convergence per iteration 
when using approximate natural orbitals here as can be 
seen in Fig. [TJ It appears that MP2 natural orbitals fol- 
lowed by those of QCISD and then MCCI all offer supe- 
rior convergence to MOs here. However neither the time 
cost per iteration nor the overhead from calculating the 
natural orbitals is taken into account. One approach is 
to consider, to three decimal places, the highest final it- 
eration energy and check how long it takes a calculation 
to first reach this energy or lower on the step after all 
states have been considered for removal. Here the final 



We now consider the carbon dimer with no frozen or- 
bitals and a bond length of R = 1.6 angstroms. Here 
the system is moving away from the equilibrium geome- 
try of R w 1.25 angstroms. We use the 6-31G* basis set 
and 200 iterations with a cut-off value of 5 x 10 -4 . In 
this system we see in Fig. [2] that the fastest improvement 
per iteration is when using MP2 NOs, then it appears 
that MCCI NOs followed by QCISD NOs improve on the 
convergence per iteration compared with MOs. However 
we find that the MP2 natural orbitals actually give the 
highest energy on the final step of —75.628 Hartree. 

We see in Table |H] that now MP2 NOs produce the 
longest time to reach the specified energy but they still 
use fewer DETs than the canonical molecular orbitals. 
MCCI NOs now perform the best, with regards to time 
and number of states to reach this energy, but there is 
not much difference between the results using MCCI NOs 
and those using QCISD NOs. 



TABLE II. Total time and number of DETs for C 2 to first 
reach E < —75.628 Hartree on a step following the consider- 
ation of all states for removal. 



Orbitals 


DETs 


Time (seconds) 


MOs 


15,045 


628 


MCCI NOs 


11,755 


377 


MP2 NOs 


13,795 


1134 


QCISD NOs 


11,953 


381 



We note that if we pick a high enough energy then MP2 
natural orbitals may give the fastest convergence and, in 
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FIG. 2. MCCI (cmin = 5 x 10 4 ) energy against number of 
iterations for C2 at R = 1.6 angstroms with a 6-31G* basis 
set when using either MOs, MP2 NOs, QCISD NOs or MCCI 
NOs. 



addition, for a single calculation the stochastic nature of 
the algorithm could affect the order of the methods. So 
for a fairer comparison we will now consider potential 
curves. As this requires numerous single-point calcula- 
tion then any random improvements or deteriorations in 
the speed of the calculation should average out and rather 
than considering an arbitrary energy as the target we will 
use a convergence criterion for each single-point calcula- 
tion. In addition to possible faster calculations and fewer 
states this should enable us to see what, if any, improve- 
ment the use of natural orbitals produce in the accuracy 
of the potential curves. 



NPE = ma,x\E 



FCI 



--&„• —mini!/, 



FCI 



Ef 



(5) 



where i ranges over all M considered points. 

One possible problem with using the NPE is that two 
curves with the same maximum and minimum error will 
have the same NPE regardless of their accuracy for the 
rest of the points. We attempt to incorporate the accu- 
racy of the other points by considering the mean squared 
value of the energy difference. Here the constant c, which 
a potential may be shifted by, is chosen to minimise the 
sum 



1 M 



(6) 
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where AE t = Ef CI - E? pprox . Setting |£ = leads to 

1 % % O r) c 



1 M 



so 



M 



M 



i = MAB 



S = ^E(A£i-MA S ) 2 = cr 
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(7) 



(8) 



This suggests the variance of the difference in energies 
a\ E as a way to quantify the fit of two potential curves 
that takes into account all the considered points and that 
the curves can be shifted by a constant without changing 
their physics. To give a quantity in units of energy we 
then use the standard deviation of AE: ctab- 



IV. POTENTIAL ENERGY CURVE COMPARISON 



A. HoO 



We now use CSFs in the main MCCI calculation, but 
we still approximate the MCCI natural orbitals using 
a DET MCCI calculation. We introduce a convergence 
check in that the calculation for each point is run until 
the maximum difference in the last three energies follow- 
ing steps where all states are considered for deletion is 
1CP 3 Hartree. Furthermore the MCCI method here is 
such that no new states are added on any iteration fol- 
lowing a step where all states have been considered for 
removal; previously this only occurred on the last itera- 
tion. This ensures that only the energies of wavefunctions 
where all states satisfy the cut-off requirement are being 
compared. We use twelve processors for the MCCI cal- 
culations except for the construction and diagonalization 
of the one matrix which is carried out in serial. 

We quantify the accuracy of the potential curves when 
FCI results are available using the non-parallelity error 
(NPE)±2 and a AE (see below). 

The NPE takes into account that a potential is defined 
only up to an additive constant so two curves differing 
only by a constant should have no error. This error is 



We consider the potential curve for the double hydro- 
gen dissociation of water at a bond angle of 104.5 degrees 
with a cc-pVDZ basis and one frozen core. We generate 
the FCI results using MOLPR01&2&2I and use a cut-off 
of c m i n = 1CP 3 in the MCCI calculations. 100 iterations 
are used to produce the MCCI natural orbitals here. 

In Fig. |3] we see that MCCI with MOs is very close to 
the FCI curve while when using QCISD NOs in MCCI 
there is a seemingly anomalous point at R = 4 Bohr. 
This may be linked to the most occupied QCISD nat- 
ural orbital having an occupation of greater than two 
(2.27 here) suggesting that the response one matrix is 
not a good approximation to the actual one matrix. Non 
physical occupations of the natural orbitals of the re- 
sponse one matrix of single reference methods has been 
suggested as a test for when multireference methods are 
required in Ref. [22- Interestingly, to four decimal places 
the largest occupancy is physical for larger R when using 
QCISD and we can see that the potential curve is again 
very close to that of the FCI. The same feature is present 
when using MP2 natural orbitals. 



FCI 

MCCI using MOs 
MCCI using QCISD NOs 




R (Bohr) 

FIG. 3. Energy (Hartree) against OH bond length R (Bohr) 
for water in a cc-pVDZ basis with one frozen core using FCI, 



MCCI (Cn 



10" J ) with either MOs or QCISD NOs. 



Quantifying the accuracy of when using these NOs for 
the the whole curve would not be useful so we instead 
first display the results for the fifteen points with R < 4 
Bohr. We include the time necessary for the calculation 
of the natural orbitals. For all the results, the time for 
the recalculation of the integrals when using MCCI NOs 
and the QCISD and MP2 calculation time are all a very 
small fraction of the total time (less than a second in 
this case) and are approximately included by using the 
time for one appropriate geometry (R = 2 Bohr for the 
first fifteen points) multiplied by the number of points 
considered. 



TABLE III. Upper part considering 
(R < 4 Bohr), lower part considerin. 
water double hydrogen dissociation in 
and ctae in kcal/mol. 



the first fifteen points 
g all twenty points for 
a cc-pVDZ basis. NPE 



Orbitals 



NPE o&e Mean CSFs Time (s) 



MOs 

MCCI NOs 
MP2 NOs 
QCISD NOs 



3.22 0.94 

4.18 1.23 

1.58 0.38 

1.35 0.32 



1507 1200 

1190 1491 

1124 795 

1050 694 



MOs 3.22 0.92 1599 1856 

MCCI NOs 4.49 1.14 1147 2263 

QCISD NOs/MOs 2.70 0.67 1256 1350 

QCISD NOs/MCCI NOs 1.35 0.37 1042 1466 



here and that the 100 iteration Slater determinant run 
does not produce natural orbitals, apart from the largest 
five, with occupations greater than 0.1 until R = 3.2 
Bohr. Furthermore the calculation to find the MCCI 
NOs has a wavefunction consisting of only a single Slater 
determinant for the smallest two bond lengths. We note 
that with only 50 iterations it was even less likely to 
produce a MCCI wavefunction consisting of more than 
one DET which is why we used 100 iterations for this 
system. In this case it would appear that the NPE is 
less good as we may not do much better compared with 
using MOs for short bond lengths and may even do worse 
yet we require more time to calculate the NOs. However 
we perhaps do better at larger bond lengths than when 
using MOs. This possible imbalance may contribute to a 
larger NPE. 

We now consider all of the twenty FCI points and use 
either MOs or MCCI NOs for the last five and do not 
consider MP2 NOs due to the slightly superior perfor- 
mance of QCISD NOs over the first 15 points. We plot 
the difference between the FCI result and that of MCCI 
using either molecular or natural orbitals in Fig. [31 There 
we see that the natural orbitals do give similar results to 
the molecular orbitals when bond lengths are small, but 
are more accurate at larger bond lengths. This shows 
how using QCISD NOs until they are unphysical then 
proceeding with MCCI NOs results in an accurate curve 
as the error has a much smaller range than with the other 
approaches. 
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n — MOs 
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R (Bohr) 



FIG. 4. Energy error (Hartree) against OH bond length R 
(Bohr) for water in a cc-pVDZ basis with one frozen core 
using MCCI (c min = 10~ 3 ) with MOs, QCISD/MCCI NOs or 
MCCI NOs. 



The upper part of Table HIl shows that QCISD NOs 
perform the best over the first fifteen points in terms 
of accuracy, time and number of CSFs followed by MP2 
NOs. It appears that, for bond lengths shorter than 4 
Bohr, the MCCI NOs perform less well with only the 
size of the final state smaller than the standard MOs and 
this is accompanied by a reduced accuracy and longer 
calculation times. We suggest that this is because the 
system is essentially well described by a single reference 



These observations are quantified in the lower part of 
Table lTiTl where it seems to be the case that MCCI NOs do 
better at longer bond lengths as the NPE is even higher 
using just MCCI NOs. The highest accuracy is achieved 
when using QCISD NOs for points with R < 4 Bohr 
then MCCI NOs for larger R where now the NPE is 1.35 
kcal/mol, half the size of when using the same procedure 
but with MOs instead of MCCI NOs. The time required 
is slightly longer when using the MCCI NOs but repre- 



sents an increase of less than ten percent compared with 
the NPE halving. The mean number of CSFs at 1042 
is a substantial reduction of the FCI space which con- 
sists of around 8 x 10 7 Slater determinants when spatial 
symmetries are neglected. 

This suggests the approach of using QCISD natural 
orbitals until the natural orbital occupations become un- 
physical then switching to MCCI natural orbitals. This 
should mean that a good approximation to the natural 
orbitals is achieved by QCISD when the correlation is 
essentially dynamic and then by MCCI when static cor- 
relation becomes important. 

The results show that ctab and NPE have the same 
behaviour with one slight difference being that QCISD 
NOs over 15 points and QCISD NOs/ MCCI NOs over 20 
points have the same NPE to two decimal places but <jae 
increases a little for QCISD NOs/ MCCI NOs showing a 
small decrease in accuracy that is not revealed with the 
NPE. Using MCCI NOs for the last five points halves 
the NPE compared with using MOs, but the cab value 
is 0.55 of its previous value. 



1. aug-cc-pVTZ 

We now increase the basis size to aug-cc-pVTZ while 
keeping other parameters the same. Fig.[5]shows that the 
potential curves behave generally as expected and when 
using natural orbitals the energy is noticeably lower as 
dissociation is approached. There are no FCI results for 
comparison but the curves and results for the cc-pVDZ 
basis suggest that the NO method should be more accu- 
rate. 




2 2.5 3 3.5 

Bond length (Bohr) 



FIG. 5. Energy (Hartree) against bond length (Bohr) of both 
hydrogens for water with an aug-cc-pVTZ basis using MCCI 
(c min = 1(T 3 ) with either MOs or NOs. 

For the first fifteen points where the QCISD response 
one matrix is physical the calculation is also substan- 
tially faster: 0.98 hours versus 4.6 hours. Furthermore 
the wavefunctions require fewer CSFs with an average 
of 1681 CSFS when using QCISD natural orbitals com- 



pared with 5744 CSFs when using MOs for the first fif- 
teen points. When using MCCI NOs for the longer bond 
lengths and considering all points we see in Table IIVI 
that the calculation is substantially faster, although the 
improvement is not as great as that seen over the first fif- 
teen points, and the mean number of CSFs is also much 
smaller. 



TABLE IV. Results for all points for the double hydrogen 
dissociation of water in an aug-cc-pVTZ basis. 



Orbitals 



Mean CSFs Time (Hours) 



MOs 5825 

QCISD NOs/MCCI NOs 1924 



7.15 
2.77 



We note that, when neglecting symmetry, the number 
of DETs in the FCI space has increased from 8 x 10 7 to 
around 7 x 10 12 when using this larger basis, which is 
approximately 9 x 10 4 as many DETs. However MCCI 
with QCISD NOs takes 3.4 times as long and with 2.4 
times as many CSFs for the points which it can be applied 
to compared with results for the method using cc-pVDZ. 
For all the points, using MOs gave a time scaling of 13.9 
and a scaling of 3.6 for CSFs compared with the cc-pVDZ 
MCCI MO calculations. When using QCISD NOs/MCCI 
NOs the time scaling was around 6.8 and about 1.9 times 
as many CSFs were required compared with this method 
using a cc-pVDZ basis. The scalings appear promising 
when compared with the growth in the size of the FCI 
space and can hopefully be further improved. 



2. Excited state 

We briefly return to water in a cc-pVDZ basis and as 
an aside we demonstrate the use of MCCI natural or- 
bitals for an excited state. Here the other types of ap- 
proximate natural orbitals considered are not available. 
We note that the current version of MCCI calculates one 
eigenvalue with the Davidson algorithm so the diagonal- 
ization routine can become unstable when dealing with 
excited states: the program may find itself in a subset 
of the CSF space so that the previous excited state of 
interest is now the ground state for example. Hence we 
only consider one geometry: the first excited state of A\ 
symmetry for water in Ci v with R = 2 Bohr. We now 
only use fifty iterations to create the MCCI NOs as many 
DETs are found for the first excited state. Furthermore 
no orbitals are frozen, however we still employ the ap- 
proximate NOs in an MCCI CSF calculation. We see 
in Fig. [6] that when using the MCCI NOs the energy is 
initially substantially higher than with MOs but rapidly 
decreases and becomes slightly lower than that due to 
MCCI MOs at convergence. The final MCCI wavefunc- 
tion uses fewer CSFs when NOs are employed here: 2308 
versus 1755. With MOs the time to convergence was 149 
seconds while only 102 seconds were needed using MCCI 



NOs however the calculation of the NOs was more in- 
volved here so the total time when using MCCI NOs was 
around 244 seconds. It would appear that fewer itera- 
tions for the calculation of the MCCI NOs may be useful 
for reducing the total calculation time. 



MCCI with MCCI MOs 
MCCI with MCCI NOs 
FCI 




20 30 40 

Number of Iterations 



FIG. 6. Energy (Hartree) against number of iterations for 
water in a cc-pVDZ basis with no frozen cores and an OH 
bond length of R = 2 Bohr using MCCI (c min = 10" 3 ) with 
MOs or MCCI NOs compared with the FCI result. 

Future work is planned using state-averaging of the 
MCCI wavefunction for a few states to reduce instabil- 
ities in the calculation and enable better calculation of 
excited potential energy curves. There is also the pos- 
sibility of other spin states being reached in the MCCI 
DET calculation and the consideration of more than one 
eigenvalue may allow better discrimination of these spin 
states. The reasonably promising result for the use of 
NOs in the calculation of an excited state should hope- 
fully be improved upon using these approaches. 



B. N 2 

We now consider the MCCI potential energy curve for 
N2 dissociation with two frozen cores in a cc-pVDZ basis. 
Here fifty iterations are used to create the MCCI NOs. 
The fifteen FCI results are gathered from Refs. I23H251 

Similar to our findings for water the MCCI run for 
small R does not result in a state beyond that compris- 
ing the occupied MOs. This occurs for both cut-offs wc 
consider, so we continue with the use of QCISD NOs 
until they become unphysical when we switch to MCCI 
NOs. This does not occur until the last FCI point (2.225 
angstroms) in this case. We see in Table IVl that the use 
of approximate natural orbitals reduces the calculation 
time and the average number of states required. The 
accuracy is also improved by the use of natural orbitals 
here. 

With a smaller cut-off, Table I VII shows that the cal- 
culation takes longer but the speedup due to the use of 
approximate natural orbitals is of a similar factor. The 



TABLE V. N 2 results with cc-pVDZ and c min = 1CT 3 . 
and (tab in kcal/mol. 



NPE 



Orbitals 



NPE (tab Mean CSFs Time (Hours) 



MOs 6.37 1.69 2909 1.69 

QCISD NOs/MCCI NOs 5.03 1.39 2478 1.10 



improvement in accuracy is not quite such a large scaling 
as for the smaller cut-off but again the approximate NOs 
have improved calculation time and accuracy. 



TABLE VI. N 2 results with cc-pVDZ and c mi „ = 5 x 10" 
NPE and uab in kcal/mol. 



Orbitals 



NPE (tab Mean CSFs Time (Hours) 



MOs 3.98 1.06 7185 7.55 

QCISD NOs/MCCI NOs 3.49 0.87 5758 4.98 

We see in Fig. [7] that the error of the MCCI results 
when compared with FCI decreases when approximate 
natural orbitals are used and when the cut-off is lowered 
from Cmin = 10~ 3 to c m i n = 5 x 1CT 4 . The reduction due 
to the smaller cut-off is greater than that due to using 
approximate NOs. 
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FIG. 7. Energy error (Hartree) against bond length (Bohr) 
for N 2 when using MCCI with a cc-pVDZ basis, two different 
cut-off values and either MOs or approximate NOs. 



C. CO 

We use a cc-pVDZ basis set to model the dissociation 
of carbon monoxide and freeze two of the orbitals. For 
MCCI the cut-off is c m i n = 5x 1CP 4 and 50 iterations are 
used for the generation of the MCCI NOs. The FCI space 
consists of 4 x 10 9 Slater determinants when neglecting 
symmetry. We find that the QCISD NOs become un- 
physical at R = 3.2 Bohr here and if we consider the 17 
points for bond length smaller than this we see in Table 
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IVHI that accuracy, time and size of the wavefunction are 
all improved by the use of NOs. Interestingly the most 
accurate curve for the first seventeen points is due to the 
MCCI NOs in contrast to the results for water and N2. 
Now the MCCI occupied NOs are never just the occupied 
MOs. The fastest calculation and fewest CSFs on average 
both belong to the calculation using QCISD NOs. 



The energy error when compared with the available 
FCI results is depicted in Fig. [5] This reveals that the 
lowest error is achieved when using QCISD NOs however 
the error increases with bond length until it becomes sim- 
ilar to that found when using MCCI NOs. The smallest 
range of errors comes from using MCCI NOs which re- 
sults in this approach having the lowest NPE and gab- 



TABLE VII. Results for the first 17 points for CO (R < 3 
Bohr). NPE and ctab in kcal/mol. 



Orbitals 



NPE o-ae Mean CSFs Time (Hours) 



MOs 

QCISD NOs 
MCCI NOs 



3.11 

2.77 
1.58 



0.89 
0.80 
0.41 



7053 
5616 
6069 



6.76 
4.74 
5.79 



The MP2 natural orbitals become unphysical sooner 
than those of QCISD for this system: the largest occu- 
pation is around 2.03 at 3 Bohr but some small negative 
higher occupations occur at shorter bond lengths. The 
MCCI point at 3 Bohr is then of poor accuracy. If we ex- 
clude this and compare over the first 16 points we have 
a NPE of 8.28 kcal/mol compared with 2.75 kcal/mol 
when using the QCISD natural orbitals. This poor per- 
formance seems to be due to the occurrence of negative, 
although small, natural orbital occupations. 

We see in Fig. [8] that the curves appear to be well be- 
haved and are close to the FCI points where they are 
available. We note that we were unable to calculate FCI 
points for larger bond lengths due to disk space require- 
ments. For the 19 points for which we have FCI results 
the NPE values in kcal/mol are 1.58 for MCCI NOs, 3.37 
for MOs, 3.40 for QCISD NOs/MCCI NOs while the 
<7A£ values in kcal/mol are respectively 0.39, 1.01 and 
0.91. It is interesting that the order of MOs and QCISD 
NOs/MCCI NOs with regards to accuracy changes in this 
case depending on whether it is quantified using the NPE 
or (Tab, although the differences are small. 
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FIG. 9. Energy error (Hartree) against bond length R (Bohr) 
for CO with a cc-pVDZ basis for MCCI (c min = 5 x 10~ 4 ) 
with either MOs, QCISD/MCCI NOs or MCCI NOs. 



We see in Table rvTTTI that . when all points forming the 
curve are considered, the use of approximate natural or- 
bitals accelerates the calculation, uses fewer CSFs and 
the potential energy curves suggest that there should be 
an improvement in accuracy. The fastest was a combi- 
nation of QCISD and MCCI NOs but the results for the 
first 19 points suggest that the most accurate results in 
this case would perhaps be due to MCCI NOs. 
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FIG. 8. Energy (Hartree) against bond length R (Bohr) for 
CO with a cc-pVDZ basis using FCI, MCCI (c min = 5 x 10" 4 ) 
with either MOs or QCISD NOs. 



TABLE VIII. Results considering all 26 points for CO. 



Orbitals 



Mean CSFs Time (Hours) 



MOs 8,208 

QCISD NOs /MCCI NOs 6,993 

MCCI NOs 7,289 



20.72 
17.88 
18.93 



We note that the 17 FCI points up to and including 
R — 3 Bohr required around 709 processor hours which 
we could approximately equate to 59 Hours when running 
on 12 processors. This is still ten times slower than the 
MCCI calculation using MCCI NOs, over the same num- 
ber of points and furthermore storage space issues meant 
that the FCI calculation could not be run to convergence 
for R > 3.6. 



V. SECOND-ORDER PERTURBATION THEORY 

The second-order perturbation scheme for configura- 
tion interaction in Ref. 15 considers an energy lowering 



a Ek =y: 



\(I\H\K)\* 
E K -{I\H\I) 



(9) 



Here \K) is the current CI wavefunction while the sum 
is over all \I) which are formed by single and double 
substitutions from \K). If a contribution from any \I) 
is greater than a threshold then these \I) are added to 
the reference space and a new wavefunction found by 
diagonalising the Hamiltonian. The process is continued 
until no new states are added to the CI wavefunction 
and then the final AEk gives an estimate of the energy 
lowering due to the neglected states. We use this scheme 
with the final wavefunction from a MCCI calculation to 
attempt to account for more of the dynamic correlation 
(MCCIPT2). For this we use an MCCI version where 
states are again added on a step following one where all 
states have been considered for removal. We note that 
the program is run for 200 iterations on eight processors 
without a convergence check here. 

If we write H = Hq + H' and have Hq \^> mcci) = 
Emcci Y&mcci)- Then for Slater determinants we have 
(I\ H \K) = (I\ H' \K), but for the non-orthogonal CSFs 
used in MCCI we need to use (I\H\K) - E K (I\ K) = 
(I\ H' \K) in the numerator to give. 



AE 



K 



E 



(I\H\K)-E K (I\K)f 
E K -(I\H\I) 



(10) 



Here all states are normalised. We use c m ; n as the thresh- 
old to consider if a state, in the PT2 scheme, should be 
added to the MCCI wavefunction. We note that for CSFs 
we use the same procedure^ of a random walk through 
the branching diagram as in MCCI. This followed by the 
removal of duplicates ensures that the CSFs are linearly 
independent but may mean that it is conceivable that 
some CSFs are neglected. 

The slowest step in the original PT2 method was check- 
ing if a prospective state was a duplicate in the set of all 
single and double substitutions or if it should be added 
to themJ^ Given the size of I as Ni then this requires 
O(Nj) operations if we check each new member against 
all previous and assume that the size of the space without 
duplicates is approximately a constant fraction of the size 
of I. As Ni is expected to be very large compared with 
the states in the MCCI wavefunction, we consider sorting 
the list of I by alpha and beta string using the quicksort 
algorithm^! This will tend to need 0(AT 7 log(JVi)) op- 
erations followed by one pass through the sorted list of 
O(Ni) to delete repeated states. We also have to delete 
any members of if in J but this is quick as K is small 
in comparison with /. The set of I can then be split 
amongst processors to calculate AE in parallel but this 



is currently implemented only in the case of Slater deter- 
minants. We note that a small test calculation with 10 
CSFs in the final MCCI wavefunction required ten times 
longer when not using the new method of removing du- 
plicates. 

We test MCCIPT2 on N 2 in a cc-pVDZ basis with two 
frozen cores and a MCCI cut-off of c m j n = 10~ 3 . The 
MCCI calculations are carried out using eight proces- 
sors. Two hundred iterations are used for each MCCI 
calculation. 

Slater determinants did not work so efficiently here: 
new states were discovered when using MCCIPT2 and 
we found that running another MCCI calculation each 
time with the reference taken as the last MCCI wave- 
function plus the added PT2 states until no new states 
were found was necessary to achieve a smooth potential 
curve. Nevertheless the use of MCCIPT2 improved the 
accuracy from an NPE of 11.01 kcal/mol for MCCI with 
PT2 states to an NPE for 6.53 kcal/mol for MCCIPT2 
while (Tab reduced from 3.12 kcal/mol to 1.86 kcal/mol. 

When using CSFs no states were found by the PT2 
procedure with a large enough contribution to be added 
to the MCCI wavefunction. This suggests that, with re- 
gards to our requirement for adding states using PT2, the 
MCCI wavefunction is, in a sense, optimum when using 
CSFs here. We see in Fig. [TO] that the MCCIPT2 curve 
appears to be of higher accuracy as at times it is diffi- 
cult to distinguish from the FCI curve on the scale of the 
graph. The plot of differences between the MCCI and 
FCI energies (Fig. ITT]) shows that the errors are much 
smaller using MCCIPT2 and a little more balanced. The 
NPE for MCCI here was similar to previous MCCI calcu- 
lations for nitrogen at 6.18 kcal/mol and this was reduced 
to 3.42 kcal/mol when using MCCIPT2. The a AE value 
was lowered from 1.83 kcal/mol to 0.92 kcal/mol by using 
MCCIPT2. 
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FIG. 10. Energy (Hartree) against bond length (angstrom) 
for N 2 with a cc-pVDZ basis using MCCI and MCCIPT2 with 
CSFs and c min = 10" 3 , compared with FCI results.— — 



The time for a a single-point MCCI calculation on 8 
processors of 200 iterations ranged from less than one 
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FIG. 11. Energy error (Hartree) against bond length 
(angstrom) for N2 with a cc-pVDZ basis when using MCCI 
or MCCIPT2 with CSFs and c min = 1(T 3 . 



minute to around 1.3 hours as R increased here. While 
the total time including the PT2 calculation on one pro- 
cessor for this proof of concept program ranged from 
less than four minutes to almost 2 hours as R increased. 
The number of states comprising the MCCI wavefunction 
ranged from around 1000 to almost 5000 as R increased 
while the states in the PT2 energy lowering calculation 
ranged from 0.4 million to 1.6 million. The accuracy 
of MCCIPT2 at c min = 10~ 3 was better than MCCI at 
Cmin =5x 10~ 4 (see Table IVl) however the time was 
longer at around 10.7 hours but this was on 8 proces- 
sors and without a convergence check. If we consider the 
time for PT2 only (4.18 Hrs) and reasonably assume it 
would be similar if used on the MO MCCI c m j n = 10~ 3 
results of Table [V] then this suggests a time of around 
5.9 Hours which would be faster than MCCI with MOs 
at c m ; n =5x 10~ 4 but not MCCI with approximate nat- 
ural orbitals at this cut-off. The results are encouraging 
and the PT2 CSF code for MCCI has room for improve- 
ment, e.g., parallelisation and more efficient calculation 
of matrix elements. 



VI. SUMMARY 

We introduced a way to approximate natural orbitals 
in MCCI and we have seen that approximate natural or- 
bitals from an MP2, QCISD or MCCI run could reduce 
the time and number of states necessary for a single-point 
MCCI calculation when using Slater determinants. We 
introduced a measure of accuracy of a potential curve 
(cab) that takes into account that the curve can be 
shifted by a constant but, unlike the non-parallelity er- 
ror, considers all points in the curve. For the curves con- 
sidered in this paper the behavior of each measure was 
usually similar although there were occasions when the 
NPE did not change but <tae did, or that the ordering of 
accuracy using the two approaches was changed for small 



differences. 

For the potential curve for double hydrogen dissocia- 
tion of water in a cc-pVDZ basis we found that if the 
QCISD or MP2 natural orbitals became unphysical the 
accuracy of the MCCI potential curve could be severely 
impacted. The results suggested the use of QCISD natu- 
ral orbitals until they had occupations greater than two 
or negative occupations, then switching to MCCI natural 
orbitals (QCISD NOs/MCCI NOs) offered the largest im- 
provement in accuracy and number of CSFs and took less 
time than when using molecular orbitals. Similar results 
were seen for the potential curve for N2 in a cc-pVDZ ba- 
sis. We noted that the MCCI natural orbitals could be 
unsuitable at bond lengths where a single reference was a 
good approximation as here the only occupied MCCI nat- 
ural orbitals were the occupied molecular orbitals when 
using short MCCI calculations. We used the approach of 
QCISD NOs/MCCI NOs for a potential curve of water 
in an aug-cc-pVTZ basis and saw good improvements in 
calculation time and the number of CSFs required. The 
scaling in calculation time compared with the cc-pVDZ 
basis was very much smaller than the increase in the size 
of the FCI space. 

For the potential curve for the dissociation of carbon 
monoxide the MCCI potential curve was most accurate 
when using MCCI natural orbitals for the points that 
we had FCI results for. The calculation time for the en- 
tire curve was a little longer than when using QCISD 
then MCCI natural orbitals but was still better than us- 
ing molecular orbitals. We note that the use of approx- 
imate natural orbitals here did not always improve con- 
vergence or reduce the error. However by using QCISD 
NOs/MCCI NOs in MCCI calculation speed and accu- 
racy were seen to be increased when compared with re- 
sults using MOs. This small sample of molecules seems 
to suggest that the MCCI natural orbitals should be 
used unless there are many MCCI natural orbitals with 
zero occupation at the start of the curve then QCISD 
NOs/MCCI NOs should be employed. 

We saw that an adaptation of a second-order perturba- 
tion scheme 1 ^ combined with MCCI (MCCIPT2) could 
run faster when using a new method to remove dupli- 
cates in the space of single and double substitutions of 
the reference. We found that at the same level of cut-off, 
the MCCIPT2 calculation with Slater determinants was 
much less efficient than that with CSFs. MCCIPT2 gave 
results with higher accuracy than the MCCI calculation 
alone for the potential curve of the dissociation of the 
nitrogen molecule. 



ACKNOWLEDGMENTS 

We thank the European Research Council (ERC) for 
funding under the European Union's Seventh Framework 
Programme (FP7/2007-2013)/ERC Grant No. 258990. 



l J. C. Greer, J. Chem. Phys., 103, 1821 (1995). 



11 



2 L. Tong, M. Nolan, T. Cheng, and J. C. Greer, Comp. Phys. 

Comm., 131, 142 (2000). 
3 J. C. Greer, J. Chem. Phys., 103, 7996 (1995). 
4 J. A. Larsson, L. Tong, T. Cheng, M. Nolan, and J. C. Greer, 

J. Chem. Phys., 114, 15 (2001). 
5 W. Gyorffy, R. J. Bartlett, and J. C. Greer, J. Chem. Phys., 

129, 064103 (2008). 
6 J. P. Coe, D. J. Taylor, and M. J. Paterson, J. Chem. Phys., 

137, 194111 (2012). 
7 J. P. Coe, D. J. Taylor, and M. J. Paterson, J. Comput. Chem., 

34, 1083 (2013). 
8 P.-0. Lowdin, Phys. Rev., 97, 1474 (1955). 
9 P.-0. Lowdin and H. Shull, Phys. Rev., 101, 1730 (1956). 
10 L. Bytautas, J. Ivanic, and K. Ruedenberg, J. Chem. Phys, 

119, 8217 (2003). 
"M. L. Abrams and C. D. Sherrill, Chem. Phys. Lett., 395, 227 

(2004). 
12 Z. Lu and S. Matsika, Journal of Chemical Theory and Compu- 
tation, 8, 509 (2012). 
13 J. A. Pople, M. Head-Gordon, and K. Raghavachari, J. Chem. 

Phys., 87, 5968 (1987). 
14 C. M0ller and M. S. Plesset, Phys. Rev., 46, 618 (1934). 
15 R. J. Harrison, J. Chem. Phys., 94, 5021 (1991). 
16 H.-J. Werner, P. J. Knowles, G Knizia, F. R. Manby, M. Schiitz, 
P. Celani, T. Korona, R. Lindh, A. Mitrushenkov, G Rauhut, 
K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhards- 
son, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, 



F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hre- 
nar, G. Jansen, C. Koppl, Y. Liu, A. W. Lloyd, R. A. Mata, 
A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nick- 
lass, D. P. O'Neill, P. Palmieri, K. Pfliiger, R. Pitzer, M. Reiher, 
T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, 
M. Wang, and A. Wolf, "Molpro, version 2010.1, a package of 
ab initio programs," (2010), see http://www.molpro.net 

17 G. D. Purvis III and R. J. Bartlett, J. Chem. Phys., 76, 1910 
(1982). 

18 J. S. Muenter, J. Mol. Spectrosc, 55, 490 (1975). 

19 X. Li and J. Paldus, J. Chem Phys., 103, 1024 (1995). 

20 P. J. Knowles and N. C. Handy, Chem. Phys. Lett., Ill, 315 
(1984). 

21 P. J. Knowles and N. C. Handy, Comp. Phys. Comm., 54, 75 
(1989). 

22 M. S. Gordon, M. W. Schmidt, G. M. Chaban, K. R. Glaesemann, 
W. J. Stevens, and C. Gonzalez, J. Chem Phys., 110, 4199 
(1999). 

23 H. Larsen, J. Olsen, P. j0rgensen, and O. Christiansen, J. Chem. 
Phys., 113, 6677 (2000). 

24 S. R. Gwaltney, E. F. C. Byrd, T. V. Voorhis, and M. Head- 
Gordon, Chem. Phys. Lett., 353, 359 (2002). 

25 G K.-L. Chan, M. Kallay, and J. Gauss, J. Chem. Phys., 121, 
6110 (2004). 

26 J. C. Greer, J. Comp. Phys., 146, 181 (1998). 

27 C. A. R. Hoare, The Computer Journal, 5, 10 (1962). 



